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Abstract 

We introduce two improvements in the numerical scheme to simulate collision and slow 
shearing of irregular particles. First, we propose an alternative approach based on sim- 
ple relations to compute the frictional contact forces. The approach improves efficiency 
and accuracy of the Discrete Element Method (DEM) when modeling the dynamics of the 
granular packing. We determine the proper upper limit for the integration step in the stan- 
dard numerical scheme using a wide range of material parameters. To this end, we study 
the kinetic energy decay in a stress controlled test between two particles. Second, we show 
that the usual way of defining the contact plane between two polygonal particles is, in gen- 
eral, not unique which leads to discontinuities in the direction of the contact plane while 
particles move. To solve this drawback, we introduce an accurate definition for the contact 
plane based on the shape of the overlap area between touching particles, which evolves 
continuously in time. 

Key words: Granular media, Discrete Element Method, Slow deformation, Contact forces. 
PACS: 83.10.-y, 45.70. -n, 45.70.Cc, 02.60. -x 



1 Motivation and model 



To model the dynamics of granular media a commonly used approach is the well- 
known Discrete Element Method (DEM) [1-3]. When applied to slow shearing [4- 
6] the computation of frictional and therefore contact forces between particles may 
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Fig. 1. (a) Sketch of the system of 256 particles (green particles) under shearing of top and 
bottom boundaries (blue particles). Horizontally periodic boundary conditions are consid- 
ered and a constant low shear rate is chosen, (b) Illustration of two overlapping particles, 
where the overlap region A between particles fully characterizes the contact force F°. 

introduce numerical errors that can be larger than the precision of the integration 
scheme used in DEM [7] . Moreover, if irregular particles are considered - typically 
polygons - the definition of the contact plane between touching particles is not 
straightforward. 

In this paper we propose an approach to compute the tangential contact force with 
the same numerical precision as the normal force. Our improved procedure is spe- 
cially suited for the case of slow shearing when large integration steps are needed 
for efficient computation, when studying e.g. the occurrence of avalanches [10] 
and the emergence of ratcheting in cyclic loading [11]. In addition, we discuss the 
drawbacks of the common procedures for computing the contact plane between 
touching particles. We will show that for polygonal (anisotropic) particles it varies 
discontinuously in time perturbing the proper convergence to the stationary state in 
stress controlled tests. To overcome this shortcoming we introduce a proper defini- 
tion which does not depend on the shape of the touching particles and is suited for 
both regular and irregular particles. 

Our model considers a two-dimensional system, as sketched in Fig. la. Each parti- 
cle has two linear and one rotational degree of freedom, being contained between 
two boundary plates shearing against each other [8, 10, 19]. The volumetric strain 
is suppressed, i.e. the position of the walls is fixed and there is no dilation. Peri- 
odic boundary conditions are imposed in the horizontal direction. The evolution 
of the system is given by the integration of Newton's equations of motion, where 
the resulting forces and momenta acting on each particle are given by the sum of 
all contact forces and torques applied on that particle, respectively. The boundary 
particles move with a fixed shear rate. 
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The integration of Newton's equation of motion is usually done with a predictor- 
corrector scheme [1, 12], which consists of three main stages, namely prediction, 
evaluation and correction. In the prediction stage one extracts for each particle the 
predicted position and acceleration of the center of mass. This is done by means of 
Taylor expansions of the linear and angular positions, yielding the corresponding 
velocities and higher-order time derivatives, as functions of the current values [12, 
14]. During the evaluation stage, the predicted coordinate of each particle is used 
to determine the contact force F^ +At at time t + At. Since the method is not exact, 
there is a difference between the acceleration f(t + At) = F£ +At /m and the value 
obtained in the prediction stage, namely Af = f(t + At) — r^(t + At). Finally, the 
difference Af is used in the corrector step to correct the predicted position and its 
time derivatives, using proper weights for each time derivative [12]. The weights 
depend upon the order of the algorithm and the differential equation being solved. 
These corrected values are used for the next integration step t + At. 

In our simulations we integrate equations of the form f — f(f,f), using a fifth 
order predictor-corrector algorithm that has a numerical error proportional to (At) 6 
for each integration step [12]. 

Typically, the particles can neither break nor deform, i.e. fragmentation is ne- 
glected. The deformation is usually modeled by letting particles overlap [1,15], 
as illustrated in Fig. lb. The overlap between each pair of particles is considered 
to fully characterize the contact: the normal contact force is assumed to be propor- 
tional to the overlap area [19] and its direction perpendicular to the contact plane, 
which is usually defined by the intersection points between the boundaries of the 
two particles [19]. 

The contact forces, F c , are decomposed into their elastic and viscous contributions, 
F e and F v respectively. The elastic part of the contact force is simply given by the 
sum of the normal and the tangential components, with respect to the contact plane, 
namely 

F e = F e n n c + Fft\ (1) 
where the normal component reads 

F e n = -kJL/h, (2) 

with k n the normal stiffness, A the overlap area and l c the characteristic length 
of the contact, which, for two particles i and j, is given by l c = r-j + tj with 

Ti = J ' Ai/2n and A,- t the area of the particle i, and similarly for particle j. The 
tangential force is considered to be proportional to the elastic elongation £ of an 
imaginary spring - also called Cundall-Stack spring [15] - at the contact, namely 

F e t = -ktt, (3) 
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where k t is the tangential stiffness and £ is the elastic elongation updated as 



£(t + At)=£(i) + ^Ai, (4) 

where At is the time step of the DEM simulation, and x7jr is the tangential compo- 
nent of the relative velocity v c at the contact point. 

Having both components of the elastic force, the tangential elastic elongation £ is 
updated according to the Coulomb limit condition \F*\ = fJ,F™, with /x the inter- 
particle friction coefficient. If the Coulomb condition is reached, particles are in the 
inelastic regime, and sliding is enforced by keeping constant the tangential force 
F\. Here the tangential elongation £ takes its maximum value ±jik n A/ (k t l c )- Oth- 
erwise, if the contact is in the elastic regime < /xF") the elongation £ can 
increase in time, following Eq. (4). 

The viscous force F v takes into account the dissipation at the contact which is 
important to maintain the numerical stability of the method. This force is calculated 
as [15] 

F v = -m r i/if, (5) 

where m r is the reduced mass of the two touching particles and v is the damping 
coefficient. 

The suitable parameters for using this model are the interparticle friction /x, the 
normal stiffness k n , the ratio k t /k n between tangential and normal stiffnesses and 
the ratio e t /e n between tangential and normal restitution coefficients. The restitu- 
tion coefficients are defined from the contact stiffness and damping coefficient [16] 
as [17] 

e n = exp {--K-q/uj) = exp ( = = ) (6) 

yjAm r k n /v% - 



for the normal component, where u> = Juj^ — rf is the frequency of the damped 



oscillator, with lu = Jk n /m r the frequency of the elastic oscillator, m r the re- 
duced mass, and 77 = v n j (2m r ) is the effective viscosity. 

We start in Sec. 2 by describing the dependence of the above numerical procedure 
on the integration step, showing that the computation of the frictional forces yields 
a numerical error larger than the rest of the calculation. To overcome this shortcom- 
ing we then describe in Sec. 3 a geometrical improvement whose associated error 
is, for any studied case, of the same order as the error of the predictor-corrector 
scheme. In Sec. 4 we address the particular case of irregular particles, discussing 
the most common ways of computing the contact plane between touching particles 
and proposing a more suitable definition. Discussions and conclusions are given in 
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Sec. 5. 



2 Choosing a proper integration step 

The entire algorithm above relies on a proper choice of the integration step At, 
which should neither be too large to avoid divergence of the integration nor too 
small avoiding unreasonably long computational time. In this section we explore 
the numerical error introduced by the frictional force in Eq. (4) and determine 
the range of proper integration steps that guarantee convergence of the numerical 
scheme. 

A typical criterion to chose the integration step At is to take a value such that 
At < t c /5 [17, 18], where t c is the characteristic duration of a contact defined 



While in several cases, one uses an integration step much smaller than the thresh- 
olds above [5], for low shear rate, very small integration steps imply a high com- 
putational effort, and so At should be chosen close to the threshold t c /5. As shown 
below, the upper threshold (fraction of t c ) below which the numerical scheme con- 
verges, strongly depends on (i) the accuracy of the approach used to calculate fric- 
tional forces between particles, (ii) on the corresponding duration of the contact 
and (iii) on the number of degrees of freedom. 

In a previous work [7] we have shown that considering values of At close to the 
upper limit t c /5 yields relaxation times depending on At for the kinetic energy 
Ek(i) = \ {jriirf + IiUjfj of each particle %, when particles are subject to a non- 
zero friction. In particular, in the absence of friction AEk attains a stationary value, 
while when friction is present it changes monotonically. 

To obtain the proper threshold as function of the parameters of our model, we con- 
sider the simple situation of two circular particles in contact, as sketched in Fig. 2, 
and study the kinetic energy of one of them under external forcing. We start with 
two touching discs, say i and j, where i remains fixed and j is subject to a force F 
perpendicular to its surface (no external torque is induced) along the a>axis. As a re- 
sult of this external force, the disc j undergoes translation and rotation. The contact 
forces are obtained from the corresponding springs that are computed as described 
in Sec. 1 and act against the external force. This results in an oscillation of disc j 
till relaxation (dashed circle in Fig. 2) with a final center of mass displacement of 
AR and a rotation around the center of mass of Ad. 



by [5,13,17], 
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Fig. 2. Sketch of the stress controlled test of two particles (discs). The particle located 
at Ri remains fixed, while the particle at Rj is initially touching particle i. The vector Rij 
connecting the center of mass of particles i and j is initially oriented 45° with respect to the 
x-axis. After applying the constant force F to disc j, the system relaxes till it reaches a new 
position (dashed circumference). Between its initial and final position particle j undergoes 
a displacement AR and a rotation A9 (see text). 

Since F is kept constant, the procedure is stress controlled. Plotting the kinetic 
energy as a function of time, yields an exponential decay as 

^.(t)=4 0) exp(-^ y ), (8) 

where t R is a relaxation time whose value clearly depends on the integration step 
At. Typically [7] the relaxation time t R depends on At for integration steps larger 
than a given threshold T t t c . We define this threshold as the proper upper limit for 
the suitable integration steps. In other words, we want to determine the upper bound 

At<T t (^k t /k n ,e t /e n )t c , (9) 

where T t (/i, k t /k n , e t /e n ) is a specific function that will be determined below. 

We start by studying the influence of the stiffness, fixing /x = 500 and e t /e n = 
1.0053. Figure 3a shows the relaxation time t R of the kinetic energy of the two- 
particle system for k n = 1, 50, 200, 10 4 and 10 8 N/m. We see that decreasing At, 
the relaxation time t R increases till it attains a maximum. This is due to the com- 
putation of the tangential force using the Cundall- spring scheme in Eq. (4) that 
yields smaller values for smaller integration steps. The stabilization of tu occurs 
when At is small compared to the natural period 1/ojq of the system. Thus, T t is 
the largest value of At for which we have this maximal relaxation time, in this case 
Tt = 10~ 3 . Further, all the curves in 3a can be collapsed by using the normalized 
integration step At/t c , as shown in Fig. 3b. From the inset of Fig. 3b we also see 
that the contact time scales as t c ~ k" 1 ^ 2 , which comes from Eq. (7). 

In the case where rotation is neglected, we obtain a constant T t = 10~ 4 as shown 
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Fig. 3. The relaxation time t r (in units of t c ) as a function of (a) the integration step At and 
(b) the normalized integration step At/t c , where the contact time t c is defined in Eq. (7). 
Here the friction coefficient is kept fixed fj, = 500 and different stiffnesses k n (in units of 
N/m) are considered. The quotient At/t c collapses all the curves for different k n . We find 

— 1/2 

t c ~ k n ' as illustrated in the inset (see Eq. (7)). As a final result one finds a constant 
T t = 10~ 3 . For other values of the friction coefficient one observes similar results. The 
relaxation time is also plotted as a function of the normalized integration step At/t c , when 
rotation is suppressed, (c) \i = 500 and different values of k n and for (d) k n = 4 x 10 8 and 
different values of [i. The dashed horizontal line in (b) indicates the relaxation time of the 
kinetic energy in the absence of friction (see text). 

in Fig. 3c. This T t value is one order of magnitude smaller than the previous one in 
Eq. (9) and can be explained by considering an 'infinite' friction coefficient in the 
rotational degree of freedom [7]. Indeed, by comparing Fig. 3c with Fig. 3a, we see 
that the relaxation time is smaller when rotation is suppressed. However, while 
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Fig. 4. The relaxation time tu (in units of t c ) as a function of the friction coefficient /x (a) 
when rotation is suppressed and (b) when rotation is considered. Here, k n = 4 x 10 8 N/m 
which corresponds to a contact time t c = 9.8 x 10~ 5 s. The normalized integration step 
At/t c = 10" 5 . 
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Fig. 3c clearly shows that tn does not depend on the stiffness k n , the same is not 
true for the friction coefficient fi, as shown in Fig. 3d. 

In Fig. 4a we observe that, in the absence of rotation, there is a change of the 
relaxation time around [i = 1, which is not observed when rotation is considered 
(Fig. 4b). This transition occurs since for larger values of /x > 1, one has F t > F n , 
and therefore the frictional force increases and drives the system to relax faster. 
Further, with rotation one observes a larger t R because there is an additional degree 
of freedom, that also relaxes. 

We also check the dependence of the relaxation time on the two other parameters 
k t /k n and e t /e n . Figure 5a shows the relaxation time as a function of the normalized 
integration step At/t c , using different stiffness ratios and with fixed /i = 500 and 
e t /e n . We point out the following. First the convergence to a stationary value of 
t R is faster for large stiffness ratio as k t /k n —> 1, yielding larger values of T t . 
Second, the stationary value of tR decreases with the stiffness ratios. This decrease 
is logarithmic, as shown in Fig. 5b. 

The dependence of tR on the stiffness ratio can be explained by taking into account 
that the larger the ratio k t /k n the larger F t . Since the tangential force F t controls the 
convergence of the numerical method, the larger k t /k n the faster the dissipation and 




c 



Fig. 5. Dependence of the relaxation time on the ration k t /k n of the tangential and normal 
stiffnesses, (a) the convergence towards the relaxation time for decreasing integration steps 
using different stiffness ration, (b) the approximate logarithmic dependence of tR on k t / k n 
(see text). Here fi = 500 and et/e n = 1.0053. No rotational constraint is here considered. 
Dependence of the relaxation time on the ration e t / e n of the tangential and normal restitu- 
tion coefficients, (a) the convergence towards the relaxation time for decreasing integration 
steps using different restitution ratios, (b) the relaxation time tR as a function of kt / k n (see 
text). Here fi = 500 and k t /k n = 1/3. No rotational constraint is here considered. 
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thus, the smaller the stationary value of Ir. The logarithmic dependence observed 
in Fig. 5b can be explained from direct inspection of Eq. (6). 

Figure 5c, showing the relaxation time as a function of the ratio e t /e n for fixed 
li = 500 and k t /k n = 1/3 can be explained similarly as in Fig. 5b. Here the 
dependence of the stationary value of t R on e t /e n is approximately exponential. 

When the Coulomb condition is fulfilled (inelastic regime), the dependence on At 
observed in all figures above, does not occur, since the strength of the tangential 
force is given by F t = jiF n . All the previous results are taken within the elastic 
regime. This indicates that the improvements in the algorithm should be imple- 
mented when computing the elastic component of the tangential contact force, in 
Eq. (4), as explained in the next Section. 



3 Improving the integration of the contact force 

Using the Cundall's spring [4], the relaxation time of the two particles only con- 
verges when At is a small fraction T t of the contact time t c . Such dependence 
results from Eq. (4),that includes At in the computation of the tangential force in 
our predictor scheme which has an error of (At) 2 . 

To overcome this shortcoming, we propose a different expression to compute the 
elastic tangential elongation £. It is based on geometric relations and does not use 
At. Our expression for £ contains only the quantities computed in the predictor 
step, guaranteeing a precision for £ of the same order as the one of the predictor- 
corrector scheme. We illustrate our approach with the simple system of two discs, as 
sketched in Fig. 2, and discuss further the more realistic case of polygonal particles. 

In Fig. 2, the elastic elongation of the tangential spring results from the superposi- 
tion of both translation and rotational degrees of freedom, £j = £^ + ^ rot \ 

For the translational contribution, the elastic elongation £ depends only on the rela- 
tive position of the two particles. In this case we substitute Eq. (4) by the expression 

ef)(t + At) = ef \ t ) + -^j-(i%(t + At) - 4(t)) ■ t c , (io) 

di + Clj J J 

where ai and aj are the radii of the discs i and j respectively, Rij is the vector 
joining both centers of mass and points in the direction i — > j (see Fig. 2). Index p 
indicates quantities derived from the coordinates computed at the predictor step. 

Considering only the rotational contribution (R\j constant), the elastic elongation £ 
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Fig. 6. (a) Comparison of the relaxation time tR (in units of t c ) when using the standard in- 
tegration scheme (squares) and the proposed improved scheme (circles). Here, we perform 
the stress control test between two particles with irregular polygonal shape, as illustrated 
in Fig. lb and rotation is neglected (see text), (b) The relaxation time tjt (in units of t c ) 
using Eqs. (10) and (11) between two discs, as illustrated in Fig. 2; for the three cases when 
considering only rotation, only translation or both, the relaxation time remains constant 
independent of the integration step (see text). 



depends only on rotation between times t and t + At: 

£ ot \t + At) = if ot \t) + (0?(t + At) - 0j(t)K, (11) 

where 0j(t) is the angle of some reference point on particle j at the predictor step 
of time t. 



Using such expression, we have shown [7] that the relaxation time t R is independent 
on the integration step. This is due to the fact that all quantities in the expressions 
above have an error of the same order of the predictor-corrector scheme. 

When considering polygonal particles, one must also take into account the shape 
of the particles. For polygons, the decomposition of the elastic elongation £ into 
translational and rotational contributions is not as trivial as for discs. This stems 
from the fact that the contact point no longer lies on the vector connecting the 
centers of mass. Thus, one should recalculate each time the position of the center 
of mass (only from translation) and the relative position of the vertices (only from 
rotation). For that, instead of using Eqs. (10) and (11), we compute the overlap 
areas between the two particles at times t and t + At and consider the dislocation 
of the corresponding geometrical centers. This yields the translational contribution 
to the spring elongation. The contribution from the particle rotation is computed by 
determining the angular change due to the rotation of the branch vector between 
times t and t + At. 



10 



Fig. 7. Illustration of the definition of contact plane (thick solid lines) between touching 
discs and touching polygons. While for discs the definition is unique, since the plane 5 C is 
perpendicular to the plain containing the center of the discs, for polygons the same is not 
true, yielding the possibility of more than one suitable choice. 

Figure 6a compares how the relaxation time varies with the normalized time step 
when the standard Cundall approach is used (squares) and when our improved ap- 
proach is introduced (circles). Clearly, the dependence on the integration step ob- 
served for the usual integration scheme disappears when our improved approach 
is introduced. Therefore, all the conclusions taken above for discs remain valid for 
polygons. In Fig. 6b we show the case of two discs for comparison. 



4 An accurate definition of the contact plane 



An additional improvement concerning polygons is the proper choice of the contact 
plane. Usually shear system consider the simple situation of spherical particles. 
Such particles yield a well defined contact plane, namely the vector perpendicular 
to the one connecting the centers of mass of the two touching particles. Figure 7 
illustrates this definition of contact plane for the case of a two-dimensional system. 
In the case of polygons the branch vectors are not parallel and therefore no unique 
definition is possible. Usually [19] one usually chooses the plane containing the 
intersection points and of the boundaries of the two touching particles. When there 
are only two intersection points there is a unique possible contact plane. However, 
due to the motion of the particles, the overlap between them may give rise to more 
than two intersecting points that yields different possible contact planes. In Fig. 7 
we illustrate a case where two different choices suit equally the condition above for 
a contact plane between two touching polygons. 
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In this Section, we propose to define the contact plane as the line passing through 
the center of mass of the overlap area having an orientation defined by an angle 
a c (with the x-axis) given by the geometrical average of the orientation angles of 
each (non-oriented) edge of the overlap polygonal with respect to the x-axis. As 
shown below such definition is unique and guarantees a continuous variation of the 
orientation of the contact plane while the polygons move, as it is the case in a model 
of irregular particles under shearing. 

More precisely, for an overlap polygon delimited by m edges, where each edge 
n = 1, . . . , m lies on a specific line y = xtana n + b n , the contact plane in the 
two-dimensional case is defined from 

y — x tan a c + b c , (12) 

where b c = y c — x c tan a c with (x c , y c ) the coordinates of the center of mass of the 
overlap polygon, and 

T m , a P 

a c = p , (13) 

with P n the length of segment n. Since the angles a n are angles between non- 
oriented lines, they can be acute or obtuse. To solve this indetermination we choose 
all the angles for the sum in Eq. (13) to be acute. Then, after computing a c we 
choose either the line in Eq. (12) or the one perpendicular to it, depending which 
one is more perpendicular to the line joining the center of mass of both particles. 

The main advantage of our definition in Eq. (13) is that while particles move the 
orientation of the contact plane varies continuously with their overlap (polygonal) 
area, as illustrated in Fig. 8 with the curve marked by bullets. On the contrary, con- 
tact planes taken from two intersection points of the two polygons present frequent 
discontinuities, as can be seen from the examples with squares in Fig. 8. Such dis- 
continuities occur since not only the size and orientation of the edges defining the 
overlap polygon vary in time, but also the number of edges. 

There is also the advantage that such a definition enables the system to behave as 
physically expected. For instance, during relaxation of a stress controlled test as 
illustrated in Fig. 9, the characteristic physical properties oscillate and converge to 
a stationary value. In Fig. 9 we plot both F t /F n and a c for both definitions. We can 
see that the convergence of the new definition is similar to a damped oscillation, 
contrary to what is observed for the conventional definition. Such deviation from 
the damped oscillation behavior occurs due to the discontinuities observed in Fig. 8, 
which act as external excitations in the damping of the particles oscillation. As a 
consequence, the relaxation time associated with the discontinuous contact plane is 
smaller than when using our definition, as illustrated in Fig. 10. 
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Fig. 8. Illustration of the orientation of the contact plane, given by angle a c as a function 
of time when two polygons overlap each other. Two different definitions of contact plane 
are considered: the usual definition (squares) and the one proposed (bullets), as given in 
Eq. (13). One sees that while the proposed definition varies continuously, the usual defini- 
tion shows discontinuities that influence the contact between touching polygons (see also 
Figs. 9 and 10). 
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Fig. 9. Time evolution of (a) the ratio Ft / F n and (b) of a c for the usual definition of the 
contact plane. The same evolution is plotted for (c)-(d) the new definition in Eq. (13). In 
each case two curves are plotted, corresponding to two different integration steps (see text). 
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Fig. 10. The relaxation time as a function of the integration step for the usual definition 
(circles) and the new definition of the contact plane (bullets). 



5 Discussion and conclusions 



We introduced a technique to improve the accuracy of the numerical scheme used 
to compute the evolution of particle systems, showing that the range of admissible 
integration steps has an upper limit significantly smaller than previously assumed. 
The new approach for computing the frictional forces is an alternative to the Cun- 
dall spring, given by Eqs. (10) and (11) and suits not only the simple situation of 
discs but also more realistic ones, where particles have polygonal shape. 

Further, our upper limit of the admissible range of integration steps was deter- 
mined for a single contact using a stress controlled test. Its dependence on the stiff- 
ness, restitution coefficient and friction coefficient was carefully addressed. When 
a larger system (N particles) is studied, multiple contacts must be taken into ac- 
count, which will yield an upper limit that we conjecture to be 1/N of the upper 
limit obtained above. 

Finally, for the case of polygonal particles, we also introduced a definition for the 
contact plane between two touching particles. Our definition is not only based on 
the intersections points of the particles in contact but also on the geometry of the 
overlap area. Thus, the contact plane varies continuously during the stress con- 
trolled tests and is unique, contrary to the usual definition. 
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